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SUMMARY 

A set of initial Cartesian coordinates, which are free of ambiguities and resonance singularities, 
is developed to study satellite mission requirements and dispersions over lon£ lifetimes. 

The method outlined herein possesses two distinct advantages over most other averaging pro- 
cedures. First, the averaging is carried out numerically using Gaussian quadratures, thus 
avoiding tedious expansions and the resulting resonances for critical inclinations, etc. 

Secondly, by using the initial rectangular Cartesian coordinates, conventional, existing 
acceleration perturbation routines can be absorbed into the program without further mod- 
ifications, thus making the method easily adaptable to the addition of new perturbation 
effects. 

The averaged nonlinear differential equations are integrated by means of a Runge Kutta 
method. A typical step size of several orbits permits rapid integration of long lifetime 
orbits in a short computing time. 

INTRODUCTION 

Several sets of averaged elements (Lorell, 1970; Broucke and Cefola, 1972; Uphoff, 1973) 
are in use for satellite lifetime studies. These usually suffer from ambiguities and resonance 
singularities for low inclinations, near circular orbit, near polar orbits, critical inclination 
resonances, and such. Moreover, it is necessary to develop the perturbation representations 
in these element coordinate systems, which often requires ingenuity and is difficult in ap- 
plication. The Cartesian coordinates have been extensively utilized and routines are avail- 
able to numerically generate most of the significant perturbations. The initial conditions 
of the Cartesian solution of the classical two-body problem have been developed for variation 
of parameters (Pines, 1961 ; Christensen, 1970;Godal and Johansen, 1968), but has not 
found wide application. With the advent of averaging as a tool for eliminating long tedious 
numerical integrations in computing solutions, this study was undertaken to reestablish the 
initial Cartesian coordinates as a useful set of parameters for orbital analysis. 
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THE INITIAL CONDITION CARTESIAN ELEMENTS 

The equations of motion of the satellite in the planetary reference frame are given by 

R = -#i — + F (1) 

where F represents the perturbation forces other than the central attraction of the principal 
body. 

The initial Cartesian coordinate parameters which describe the motion are given in terms of 
the position and velocity vectors in Cartesian coordinates by 


Rq = -fR + fR 


( 2 ) 


where f, g, g, f are given as functions of the difference in eccentric anomaly, 0, as 


f =1- 
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The relationship of the time increment to 0, in the absence of perturbations, is given by a 
form of Kepler’s equation 
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(4) 


t-t 


0 




r d 

0 -sind +— sin^ 
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(1 - cosfl) 


The differential equations for the variation of and under the action of F are generated 
in Pines (1961) from the conditions that 

d 

— R = 0 
dT 


and 



(5) 


Thus, we have 

^R„=g^R-g^R-gF 


d . . 

— R„=-fR + fR + fF 
dr ° ^ 


( 6 ) 


Following Godel and Johansen (1968), we choose for the perturbation equation of the 
difference in eccentric anomaly, 6 


6 


= 0 


(7) 


This serves to eliminate mixed secular terms (see Christensen, 1970) from the perturbation 
derivatives of f, g, f, g, that arise in Pines (1961) where the time from the initial time was 
assumed unperturbed. 




(8) 


We can replace the vectors R, R from (6) using the inverse of (2) 

R = fRo+gRo 

R = fRo+gRp 


(9) 


^ Rq = (g/ - fg,) Ro + (ggr ' 8r8) Ro ‘ 


^Ro=(f,^ff)Ro^(M-f^)Ro + fF 


(9a) 
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An important stability condition on (9a) is derived. We have 

fg-gf= 1 


It follows from the perturbation derivative of (10) that 

f^-f^ = -g/+fg. 


( 10 ) 


( 1 !) 


Thus the coefficient of in the d/dr expression is the negative of the coefficient of 
Rg in the d/dr R^, expression in Equations (9a). This will reduce the computing work in 
the averaging process. 


Since the perturbation of the transit time is not zero, we require a differential equation for 
the time. Using (3), (4), and (5), we obtain the perturbation differential equations for the 
time, 



1 

V^T 



- sin 0) + sin 0 - 
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(12) 


where 


a. =2a^ R F 


= R • F 


The perturbation derivatives of f, g, f, and g are given below: 
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(13) 
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where 


-cos0)a^- 


V^sinfl 

nAT 



(13a) 


We now proceed to the averaged differential equations. 

THE AVERAGED EQUATIONS OF MOTION 

The principle here is to replace the nonlinear differential equations for R^fT) and R^,(t) by 
their average value over a single period of the reference orbit as defined by Rq(t, ) and 

Let the period of i 'erence orbit be T(t, ), then for 


T T 


(14) 


Thus, 
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where 


r *io 

— *I-cosO + — cos 9 -I- — sm( 
* * \'ia 


(15. 


Equation (15) couk) be integrated analyticaUy, using Fourier series in 9. However, this 
would require a representation of F in d and greatly burden the introduction of additional 
perturbations. Moreover, mathematical resonances would appearand require special 
techniques for each resonant term. In this, wc follow the lead of Uphoff ( 1973) and adapt 
Gaussian quadratures s the technique few evaluating the integrals. 

We note that F is a vector in three space, and recalling ( 1 1 ), we need to evaluate only 1 0 
int^rals. The average equations become: 

=a,Ro + a,R^+C 
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(cont.) 


To obtain the averaged F, G vectors, it is necessary to specify the perturbations. For the 
purpose of this study, we consider accelerations due to tesseral and zonal harmonics for a 
rotating planet, third body forces, such as the sun and moon, and drag. The detailed 
equations are given in Appendix A. It should be noted that where F is a function of R, 

R and t, for the purposes of averaging, these functions are given by the two body equations 
and Kepler's equation given in (17a) through (I7d) referenced to the r time at the midpoint. 

We now consider the numerical integration of the averaged equations. We propose to use 
a Rui^e Kutta method with a r step equal to several periods. It must be borne in mind that 
the time variable is r and not time. For each evaluation of the seven derivatives ( 1 6), 

at a specific r time, Tj, we must carry out the averaging procedure. Using Rq(t.) R^(r.) as 
the reference orbit, we compute the 10 integrals given by (16a). 

The Runge Kutta solution will produce Rj(r), Rj,(r) and (t - 1^,) (t). To compute an 
ephemeris of the state, R (t) and R (t), we proceed as follows: 

Let the period of the R„(v), R^fr) orbit be T (t). 


where 


T(t)= — (a(T))3/^ 
P 



R„(r)R„(r) 


Let N be the integer part of t/T (t), and d is given by, 

= 2v N + a 

Using Newton’s method, solve Kepler’s equation for a 


At = t-NT = 


v5T 


(o - sin o) + r^a'^ sin o + 


d^a( I - cos a) 

y/ir 


(I7a) 


(17b) 


(17c) 
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For any position in the N + I revolution, where the incremental eccentric anomaly, 6, lies 
between 0 and 2 tt. 

Let 


t = t(T| + 


r(T) d(T)(l - cos a) 

o-sina+— sino+ (l+a^(r)) 


y/ia 


R(t)=f(«)R„(r) + g(a) R„(r) 


R(t) = f(a) R„(r) + g(tt)R„(r) 


(I7d) 
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APPENDIX A 


We consider several typical |x:rturbations. 

Tesseral and Zonal Harmonics. 

Let the central body rotate about its polar axis, k (t) with a uniform angular rate cj. If the 
inertial Cartesian coordinates of the vehicle is R, then the unit vector to the vehicle in planet 
fixed coordinates is given by, 

R 

u =k(t) — 
r 


R 

s =i(t)* — 
r 


(A.l) 


R 

w=j(0 — 
r 

where i (t) and j (t) are orthogonal reference axes fixed in the rotating body equation, 
perpendicular to the polar axis k (t), expressed in the inertial Cartesian system of the 
vehicle. 

Following Pines, ( 1973). the acceleration in the inertial system is given by 

R(t) 

F| =a,i(t)+ajj(t)+ajk(t)+a^ — — (A.2) 

For the purpose of completeness, we list the a. coefficients for J, . J , . . and . 

a,( ) = ttjO j ) “ ^ “ “il-* J ) ~ *1 1-*4 I = ^ ® 

ftal 3Jj (I - 5 U-) 




2r" 


a4(J3) = - 


fia^ 5Jj(Jtu-7u^) 


(A.3) 
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(A3) 

(toat.) 


3Jj (5u^ - I) 
“3<J3> = -— ^ 


nil 5J^ (7 u* - 3 u( 
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2r* 

fiaj6 

&,(t Sjj) = — (wSjj + s^22^ 


paJ-6 


«,(C,,.S^) = 


(sS 


22 


■wCj,) 


a,(C^,.S^) = 0 


lial 15 

a^(Cjj.S^) = -— — -(C^(s^ -w^) + 2swSj2) 

Third Body Acceleration. 

Let the gravitational mass coefficient of third body be u.. Let the ephemeris of the third 
body, R (t), be given with respect to the central planet, about which the vehicle is orbiting, 
then the indirect acceleration is given by 



(A.4) 


where 


r . - K • • 

Cl Cl 


r = IR-R,;! 

VI Cl 


(A.5) 


Atmospheric Drag. 

Let p be the air density, given as a function of the vehicle position. R, and the time. t. The 
drag acceleration is given by 


F 


3 


-jP(R,t) 


vCnS . 

R 

m 


(A.6) 


no 



where 


V = IRI 

Cp ' coefficient 

\w = vehicle mass 
s = effective vehicle arag area. 


^ QVALn 
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NONLINEAR COUNTER EXAMPLE FOR BATCH 
AND EXTENDED SEQUENTIAL ESTIMATION ALGORITHMS 


B. T. Fang 

Wolf Research and Development Corporation 
Riva-dale, Maryland 


A simple nonlinear example is presented which shows the welldcnown iterated batch least- 
squares and extended sequential estimation algorithms may converge to different estimates. 
For this example one may even say the extended sequential algorithm converges to the 
“wrong” value. 
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ORBIT DETERMINATION IN THE PRESENCE OF 
UNCERTAIN ATMOSPHERIC DRAG 


B. Tapley, D. Dowd, and B. Schutz 
University of Texas 
Austin, Texas 


Uncertainties in the knowledge of the atmospheric density in the associated drag parameter 
constitute one of the primary limitations on the accuracy on which the orbits of near earth 
satellites can be determined and predicted. In most orbit determination programs, the effect 
of uncertainties in atmospheric drag are determined by adopting a standard atmosphere and 
estimating the drag parameters, 0 C^ A -j- m. 

However, for most missions, Cp and A vary and the standard atmosphere will contain errors. 
Each of these factors will lead to errors in the orbit determination and prediction operation. 
In this presentation, an approach for estimating the drag parameter, the effective satellite 
cross sectional area in the atmospheric density simultaneously with the satellite state, is 
described. 
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